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Entropy serves as a central observable in equilibrium thermodynamics. However, many biological 
and ecological systems operate far from thermal equilibrium. Here we show that entropy produc- 
tion can characterize the behavior of such nonequilibrium systems. To this end we calculate the 
entropy production for a population model that displays nonequilibrium behavior resulting from 
cyclic competition. At a critical point the dynamics exhibits a transition from large, limit-cycle like 
oscillations to small, erratic oscillations. We show that the entropy production peaks very close to 
the critical point and tends to zero upon deviating from it. We further provide analytical methods 
for computing the entropy production which agree excellently with numerical simulations. 
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The study of complex systems with a large number 
of interacting particles requests global observables that 
characterize their behavior. Modern statistical mechan- 
ics has successfully identified, interpreted and applied 
such observables for equilibrium systems. One of these 
observables is the entropy which allows for predictions 
of a system's behavior through the second law of ther- 
modynamics - an isolated system's entropy cannot de- 
crease. Identifying similar principles for non-equilibrium 
systems, however, proves elusive. Neither a characteris- 
tic global observable nor a universal principle have been 
identified in a general way. While also in nonequilibrium 
the entropy production has been proposed as a useful 
observable [U [5], and different principles governing its 
behavior have been suggested O |4] problems arise from 
different employed definitions of entropy and approaches 
to nonequilibrium dynamics [21 [51 [S] . 

In this Letter we demonstrate that entropy produc- 
tion can successfully characterize ecological systems with 
cyclic competition. Ecological systems display a wide va- 
riety of nonlinear and nonequilibrium behavior. Random 
interactions between individuals and the finitcness of the 
population lead to intrinsic stochasticity. Nonequilib- 
rium results when interactions between individuals of dif- 
ferent species include cyclic dependencies where a species 
Ai benefits from and suppresses a species A^, while A2 
benefits from and suppresses a species A3 and so on, with 
some species A^ of the resulting chain benefitting from 
and suppressing species Ax. Such cycles can lead to er- 
ratic or limit-cycle oscillations in the steady state of the 
population dynamics [7lI16j. Experimental observations 
of cyclic dynamics and corresponding oscillations have, 
amongst others, been documented for mating behavior 
of lizards in costal California [17| and in microbial labo- 
ratory communities |18j . 

The dynamics of ecological systems can be conve- 
niently described as a Markovian stochastic process 



through a master equation. 
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in which Pi (i) denotes the probability of finding the sys- 
tem in a certain state i at time t and w^ is the transition 
probability from state j to state i. The associated mean 
entropy production S of the system follows as 
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Equation ^ can be obtained through considering the dif- 
ference between forward and backward entropy per unit 
time of the stochastic process, Equation ((I]) [19 . Equa- 
tion ((2| follows also as the temporal derivative of the sys- 
tem's Gibbs entropy together with a term describing the 
total increase of thermodynamic entropy in the reservoirs 
to which the system is coupled [20] ■ For steady states de- 
fined by dtPi{t) — 0, as we consider in this Letter, the 
entropy production simplifies to 



S = 
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It follows from Equation (pi) that the entropy produc- 
tion vanishes if and only if the system obeys detailed bal- 
ance, w^-Pj = w*Pi. Indeed, detailed balance represents 
the notion of thermodynamic equilibrium in the frame- 
work of the master equation. Cyclic population dynamics 
violates detailed balance; the computation and discussion 
of the associated entropy production is the scope of this 
Letter. 

Consider a simple model for cyclic population dynam- 
ics of three species A,B, and C. Interactions are formu- 




FIG. 1: (color online) Probability distributions (fc = 1, N = 
100). (a) For a mutation rate m — 0.003 smaller than m^ = 
k/{2N) the probability distribution is concentrated near the 
edges and particularly near the corners of the phase space, (b) 
A mutation rate m = 0.1 larger than rric leads to a gaussian 
distribution around the center. 



lated as chemical reactions: 



AB 
BC - 

CA 



■ AA, 

BB, 

>CC. 



A^^B, 
B^^C, 

C^^A. 



(4) 



The reactions on the left describe cyclic competition: A 
outperforms B but is beaten by C, and C is taken over 
by B in turn. The reactions on the right correspond to 
spontaneous mutations between the three species. 

The population model defined by the reactions Q ex- 
hibits a critical mutation rate that, in the resulting non- 
equilibrium steady state, delineates large oscillations in 
the species densities from only small ones. Let us intro- 
duce this critical mutation rate first. The reactions Q 
conserve the total number N of interacting individuals. 
The densities a, 6, and c of species A, B, C therefore sum 
up to one: a + b + c = 1 and the population's state space 
is the simplex 5*3, see Fig.[l] Numerical simulations indi- 
cate that small values of the mutation rate m lead to large 
oscillations between the densities of the three species; the 
probability distribution is highest close to the corners of 
the simplex [Fig. flla)]. Large values of m, on the con- 
trary, lead to an approximately gaussian probability dis- 
tribution around the simplex center [Fig. [l|b)]. Erratic 
oscillations occur at small amplitudes [9^. 

The system's behavior can be analytically described by 
an approximate Fokker-Planck equation. A systematic 
expansion in the system size TV yields an equation for the 
temporal evolution of the probability distribution P(s, t) 
of the densities s = {a,b) at time t: 



dtPisJ) = -d^[a,{s)Pis,t)] + -d,d,[P^j{s)P{s,t)] 



(5) 



in which the indices i, j run from 1 to 2; the summation 
convention implies summation over them. The density c 



follows as c = 1 — a — fe. The coefficients read 

ai{s) = [m(l - 3si) + fcsi(sj+i - Sj+2)] , 

P^^{s) =N^'^ [m(l + Si) + kSi{Si+i + Si+2)] , 

Pijis) =- N^^['m{si + Sj) + ksiSj] for i ^ j , (6) 

where the indices are understood as modulus 3 and 
S3 = c. The terms containing a describe the determinis- 
tic part of the temporal evolution. In the absence of fluc- 
tuations, the reactions for cyclic dominance lead to neu- 
trally stable oscillations around the internal fixed point 
s* = (1/3, 1/3, 1/3), while the spontaneous mutations 
render the internal fixed point stable. Demographic fluc- 
tuations are, for large system sizes N, inversely propor- 
tional to \/N and enter the Fokker-Planck equation ([5| 
through the terms containing /3. They induce a stochas- 
tic drift away from the internal fixed point towards the 
boundaries of the phase space. The Fokker-Planck equa- 
tion ([5]) shows that the competition between the deter- 
ministic and the stochastic effects leads, at a critical mu- 
tation rate 771^ — k/{2N), to a uniform probability distri- 
bution. Certain deviations from the uniform distribution 
occur near the phase space boundaries where the dis- 
creteness of the phase space becomes relevant and the 
continuous formulation through the Fokker-Planck equa- 
tion does not hold. For small mutation rate, m < rric, 
fluctuations dominate and drive the system towards the 
boundary. In the absence of mutations the corner states 
are absorbing and the system goes extinct [5T]. An ar- 
bitrary small mutation rate, however, leads to sustained 
species coexistence and oscillations. In the opposite case, 
when m > rric, the deterministic dynamics centers the 
probability distribution around the internal fixed point. 

The cyclic population dynamics yields a non- 
equilibrium steady state that is characterized by oscil- 
lations, large or small, around the internal fixed point. 
What is the resulting entropy production and how does 
it relate to the regimes of small, critical and large muta- 
tion rates outlined above? 

To tackle this question we have carried out extensive 
numerical simulations of the stochastic system employ- 
ing the Gillespie algorithm 22J. Throughout our simula- 
tions we have considered fc = 1 which defines the time- 
scale. Numerical results from computer simulations of 
the stochastic system show that the entropy production 
peaks at a certain value mmax of the the mutation rate 
(Fig. pi) . The value m,„ax approximately equals the crit- 
ical mutation rate, TOmax ~ itt-c (Fig. ^ inset). Small 
deviations from this behavior arise for the probability 
distribution at the critical mutation rate is not uniform 
near the boundaries as mentioned above. 

Analytical understanding of the entropy production in 
the regimes of small, critical, and large mutation rates 
is feasible through the Fokker-Planck equation ((sl). To 
this end we employ a continuous version of the entropy 
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FIG. 2: (color online) Entropy production in the steady state 
for different system sizes (Z, N = 2; M, N = 5; o, N = 10; •, 
Af = 20; A, iV = 50; , N = 100; v, A^ = 200; T, iV = 400). 
The entropy production vanishes for very high and very low 
mutation rates and exhibits a maximum at an intermediate 
value rrimax- The value rramax is near the critical mutation rate 
as shown in the inset where the black line indicates rric+O.OOl, 
and red circles represent data obtained from simulations. 



production ([3]), 
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where integration is over all states r , s of the phase space. 
The entropy production in form of Eq. ([7| can readily 
be evaluated at the critical mutation rate rric. The prob- 
ability distribution is uniform according to the Fokker- 
Planck equation ([5]); we obtain 



S^^^^ = —kN [12\n{N) - 13 



61n(4)]. (8) 



For moderate and large N the term iVln(A^) on the right- 
hand side dominates the entropy production. Stochastic 
simulations confirm this behavior (Fig. [2]). 

In the regime of large mutation rates, m > rric, we need 
to calculate the probability density in the steady state to 
compute the entropy production. We obtain the proba- 
bility density by using polar coordinates (r, 0) centered 
at the internal fixed point. We then simplify the Fokker- 
Planck equation ([5| through a van-Kampen approxima- 
tion for the coefficients (|6|: the latter are approximated 
by their values at the internal fixed point. The resulting 
Fokker-Planck equation is then solved by the gaussian 
distribution 



P{r,^) = 



27Ta' '^^P \ 9^2 
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where cr = y ^omN ■ "^^^ entropy production follows from 
Eq. (It]) where the integral is evaluated by setting the 
upper boundary of the integral to oo and an average over 



FIG. 3: (color online) Entropy production in the limiting 
cases m <^ rric (a) and m 2> rric (b). Analytical results 
(black lines) agree excellently with simulations { , N = 30; 
A, iV = 100; n, TV = 200; O, iV = 400). The data confirm 
that the entropy production is proportional to the squared 
system size A''^ for small mutation rates and proportional to 
N for large mutation rates. The simulation results further 
confirm that the entropy production decays as m for jn — )■ 
and as 1/m for m — ^ oo. 



the angular dependence is taken 
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This result agrees excellently with numerical simulations 
[FigjSJ^b)]. The entropy production (10 1 depends linearly 
on the system size N. This behavior arises because the 
typical area in phase space explored by the dynamics is 
proportional to a^ ^ \/N and thus contains TV^ct^ ^ N 
states. The continuity approximation employed in the 
Fokker-Planck equation ([5]) holds for arbitrary large ttt-, 
since the width a of the probability distribution ^ re- 
mains finite as TO -^ oo. 

Expanding ( 10 1 for large values of m results in 
k'^N/{9m). The entropy production vanishes 
Indeed, only the cyclic dynamics at 



as m/k increases, 
rate k underlies the nonequilibrium behavior and there- 
fore entropy production, while the mutations at rate m 
obey detailed balance. 

When the mutation rate is small, m <C rric, the proba- 
bility distribution is concentrated near the boundaries of 
the phase space [Fig. [I] (a)]. The dynamics occurs pre- 
dominantly along the boundary and can therefore be ap- 
proximately described by considering only the boundary 
states. Because of the threefold symmetry it suffices to 
regard only one edge of the simplex with periodic bound- 
ary conditions. The concentration x of one of the three 
species increases along this edge from to 1 such that 
the cyclic dynamics drives the system to a; = 1. The 
deterministic part of the dynamics is given by 



dfX = m{l — 2x) — k{x — x ) 



(11) 



and features a fixed point at x* = [(2m + k) — 
•\/4m^~+~A?]/(2fc). In the range of m ^ mc, this 
fixed point is closer to 1 than the distance 1/iV be- 
tween two discrete states. We conclude that fluctua- 
tions will cause the system to exhibit a constant cir- 
cular current in the steady state. The probability dis- 
tribution P{x) for X G [jj^^ ~ W^ '^^^ therefore be ob- 
tained as solution to the Fokker-Planck equation — 
—dri;{[m{l — 2x) — k{x — x'^)]P{x)] where fluctuations 
have been ignored: 



P{x)^U- 



1 



— 2mx — kx + kx^ 



(12) 



with a normalization coefficient M. To determine the 
probability Pq of a corner state, which turns out to be fi- 
nite, fluctuations have to be included. Pq can be obtained 
using the master equation and the values of P{x = 1/N) 
and P{x = 1 — 1/-/V). The normalization Af follows from 



i-i 



P{x)dx+-Po 



= 1. 



(13) 



The factor 3 arises because the phase space simplex pos- 
sesses three edges. For moderate and large system sizes 
N we obtain Pq = M /{2m) which dominates the left- 
hand side of Equation (13), such that N — Nm/3. The 
resulting probability density can again be inserted into 
(It]) to provide an analytical result for the entropy pro- 
duction in the regime m <C ^ ^ 1: 



S^, 



m<^mc. 



mN \ii{m/k)^ 



(14) 



in perfect argeement with simulations [Figp^(a)]. The 
entropy production for small mutation rates is propor- 
tional to the squared system size. Decreasing m lowers 
the entropy production in proportion because mutations 
are the process that restart the cyclic dynamics once a 
corner state has been reached. Mutations therefore limit 
the dynamics to a timescale proportional to m. 

In conclusion, we have examined the global entropy 
production in the steady state of a cyclic population 
model. At a critical mutation rate the system under- 
goes a transition from large oscillations along the phase 
space's boundary to small erratic oscillations around an 
internal fixed point. The entropy production peaks very 
near the critical mutation rate and decreases to zero away 
from it. We believe that, in a similar manner, the en- 
tropy production can yield valuable information about 
the nonequilibrium steady state of other stochastic sys- 
tems. Indeed, in a recently studied model, because of 
a non-fixed system size and the extensivity of the en- 
tropy production, the slope of the entropy production 
peaks near a critical point [25]. Within our approach of 
a fixed system size we have investigated a stochastic sys- 
tem that displays a supercritical Hopf bifurcation. We 
found that the entropy production predicts the scale of 



the critical mutation rate: it peaks near the Hopf bifur- 
cation, at a mutation rate of about 1/4 of the critical 
one [23] . Understanding the certain discrepancy between 
the maximum and the critical value may open a route to 
more general understanding of the role of entropy pro- 
duction. Because of the universality of the Hopf bifur- 
cation we conclude that our approach is valid for a wide 
class of non-equilibrium systems, namely those that ex- 
hibit a transition from small, erratic oscillations to limit- 
cycle-like ones, including systems with spatial degrees of 
freedom [24j . 
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Entropy production can characterize the behavior of a broad class of nonequihbrium systems. In this Supplemen- 
tary Material we underpin this point through consideration of a nonequilibrium stochastic system that exhibitis a 
Hopf bifurcation. We show that the entropy production peaks in the vicinity of the bifurcation, where the behavior 
changes from noisy, erratic oscillations to larger limit-cycle oscillations. We conclude that the entropy production gen- 
erally characterizes the behavior of systems with limit cycles that fall into the universality class of the Hopf bifurcation. 

Consider a stochastic system with species A, B, C and empty sites that obey the following reactions: 



AB 
BC 



CA 



A(d, 
B(Z), 
C0, 



A0 - 

B<Z) 

C0 



AA, 

■ BB, 

CC, 



A< 

B 

C 



B, 

■c, 

> A. 



(15) 



The reactions with rates k and I can represent cyclic dominance of three species [11]. The corresponding deterministic 
rate equations have first been proposed and analyzed by R. M. May and W. J. Leonard [16]. The reactions with rate 
771 describe spontaneous mutations between the three species. 

The deterministic equations for the temporal evolution of the concentrations a, 6, c of species A, B, C follow from 



the reactions ( 15 ) as 



dta — a[l{l — p) — kc] + m{b + c — 2a) , 
dtb = b[l{l - p) - ka] + m{a + c-2b), 
dtc = c[l{l - p)- kb] + m{a + b-2c). 



(16) 



Linear stability analysis reveals the existence of a reactive fixed point at (a*,5*,c*) = 1/{H -I- k) ■ (1,1,1). This 
fixed point changes its stability at a critical mutation rate TOc = kl/6/{3l + k) from a stable spiral point (above rric) 
to an unstable spiral point (below rric). Further analysis that takes the leading nonlinearities into account reveals 
that a supercritical Hopf bifurcation arises at rric. Above the critical mutation rate, for m > rric, the stochastic 
system performs noisy erratic oscillations around the reactive fixed point. The steady-state probability distribution is 




(b) 




FIG. 4: (color online) Steady-state probability distributions for k = I = 1, N — 100, projected on the plane spanned by 
(a, b, c) = (1, 0, 0), (0, 1, 0), (0, 0, 1). (a) For a mutation rate m = 0.022 smaller than rric ~ 0.042 the probability distribution is 
concentrated along the limit cycle, (b) At the critical mutation rate m^ a broad probability distribution centered around the 
reactive fixed point arises, (c) A mutation rate m — 0.062 larger than rric leads to a narrow, gaussian distribution around the 
reactive fixed point. 



■Co 

o 
O 

o 

CD 



0.15 




0.05 



0.001 0.01 0.1 

(mutation rate) m 



FIG. 5: (color online) Entropy production in the steady state. Results for different system sizes (_, A'' = 100; ■, A*' = 200; o, 
N = 400; *, N — 800) show that the entropy production peaks at a mutation rate rrimax ~ 0.01 near the critical mutation rate 
rric ~ 0.042. The entropy production vanishes both for smaller and higher mutation rates. 



approximately gaussian around the reactive fixed point, see Fig.l4](c) for a projection of the system's steady state onto 
the simplex spanned by the densities a, b, c. Below the critical mutation rate, for m < rric a stable limit cycle forms. 
The stochastic dynamics leads to noisy trajectories along the limit cycle, see Fig. l4](a). At the critical mutation rate, 
as the linear terms in the deterministic equations vanish, a relatively broad, non-gaussian probability distribution 
centered at the reactive fixed point arises [Fig.l4](b)]. This behavior is similar to the one recently reported in Ref. [13] 
where higher order nonlinearities render a spiral point stable while the linear terms vanish. 

We have performed extensive stochastic simulations of the stochastic system defined by the reactions ( 15 ). In these 
simulations we have left the rates k,l constant at fc = / = 1, defining the time-scale, and systematically varied the 
mutation rate m as well as the system size N. In principle, a divergence in the entropy production can arise when 
the system reaches the boundary of the phase space. However, because the probability of these boundary states is 
exponentially suppressed, this effect can be ignored. 

For all considered system sizes the resulting entropy production exhibits a maximum near the critical mutation 
rate ttIc, see Fig. [5] For system sizes above about N = 200 the maximum of the entropy production arises at a value 
'Timax ~ 0.01, about 1/4 of the value of the critical mutation rate rric ~ 0.042. For mutation rates much smaller and 
much larger than Wmax the entropy production tends to zero. The system's behavior therefore resembles the one 
reported in the main part of this Letter, underpinning the general usefulness of entropy production in characterizing 
nonequilibrium steady states. Understanding the certain discrepancy of rrimax and rric will yield further insight into 
the relation between entropy production and critical nonequilibrium behavior. 



